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Abstract 

Molecular dynamics simulations of ionic solutions depend sensitively on the force fields employed 
for the ions. To resolve the fine differences between ions of the same valence and roughly similar 
size and in particular to correctly describe ion-specific effects, it is clear that accurate force fields 
are necessary. In the past, optimization strategies for ionic force fields either considered single- ion 
properties (such as the solvation free energy at infinite dilution or the ion-water structure) or ion- 
pair properties (in the form of ion- ion distribution functions). In this paper we investigate strategies 
to optimize ionic force fields based on single-ion and ion-pair properties simultaneously. To that 
end, we simulate five different salt solutions, namely CsCl, KCl, Nal, KF, and Csl, at finite ion 
concentration. The force fields of these ions are systematically varied under the constraint that the 
single-ion solvation free energy matches the experimental value, which reduces the two-dimensional 
{fj, e} parameter space of the Lennard Jones interaction to a one dimensional line for each ion. From 
the finite-concentration simulations, the pair-potential is extracted and the osmotic coefficient is 
calculated, which is compared to experimental data. We find a strong dependence of the osmotic 
coefficient on the force field, which is remarkable as the single-ion solvation free energy and the 
ion- water structure remain invariant under the parameter variation. Optimization of the force field 
is achieved for the cations Cs"*" and K+, while for the anions I" and F~ the experimental osmotic 
coefficient cannot be reached. This suggests that in the long run, additional parameters might 
have to be introduced into the modeling, for example by modified mixing rules. 

PACS numbers: 
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I. INTRODUCTION 



Aqueous electrolyte solutions are of fundamental importance not only in physical chem- 
istry but also for biological function and technological applications. In biology, the presence 
of ions, specifically of the monovalent ions K"*", Na^, Cl~, has significant effects on the stabil- 
ity, structure, and function of nucleic acids and proteins and the regulation of biomolecular 
processes [l, 2, 3]. In technological applications, ions can play an important role in chem- 
ical reactions by influencing their rates, as well as in controlling the solubility of various 
co-solutes [4, J For salt concentrations larger than ~ 10 mM salt effects are typically 
ion-specific even for simple bulk properties such as the osmotic pressure, which, in turn, can 
be highly relevant for transport and function in biomolecular systems The molec- 

ular understanding and prediction of the complex and often context-dependent effects of 
aqueous electrolytes poses a challenging task to the scientific and, in particular, theoretical 
community 0, 31 . 

The successful molecular modeling of ionic effects typically involves computer simulations 
in which the ionic and water degrees of freedom are explicitly resolved and evolved by a set of 
effective, classical interactions: the simulation force field. Quite commonly used force fields 
are pair-wise additive and non-polarizable to keep the parameter space small. On that level, 
the atoms are characterized by (partial) Coulombic point charges g^, excluded- volume radii, 
and dispersion attraction strengths. In standard protocols, the non-electrostatic interaction 
for atoms i and j at a distance is modeled by a pairwise Lennard- Jones (LJ) interaction 
of the form 
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cr,; 



(1) 



with two free parameters, the interaction length and the energy scale Sij, per pair of 
atoms. The whole set {cij, e^j, g^} with i, j = 1, M, defines the total force field behind the 
nonbonded inter- and intra-molecular molecular dynamics (MD) interactions for M atomic 
species. Typically, the vast number of parameters is reduced by using heuristic mixing rules 
for the cross interactions {i j) so that the only remaining parameters are the diagonal 
coefficients an and en. The common mixing rules are Eij = ^EnEjj^ and either aij = 
{aii + ajj)/2, constituting the Lorentz-Berthelot (LB) mixing rules, or o"jj = y/dnajj, defining 
the geometric mixing rules . It is assumed that the force fields take implicitly into account 
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the polarizability, as well as many-body effects. In fact, potentials which account for electric 
polarizability do not strictly seem to be required for modeling ion pairing: It has been shown, 
that for mono- and divalent ions even the first hydration shell is not significantly polarized 
compared to water molecules in the bulk 9|], though it should be kept in mind that a force 
field with more parameters has the principal possibility to be more accurate. For water 
usually simple point charge models (e.g., TIP3P or SPC/E) are used in which oxygen and 



hydrogen atoms are resolved 
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The latter are connected by rigid intramolecular bonds 



and carry partial charges optimized in such a way that a few important water properties 
(density, structure, surface tension, dielectric constant) are well reproduced. 

Throughout the years, there have been numerous attempts to systematically optimize 
ionic force fields. In principle, to properly describe the interactions between ions and be- 
tween ions and water, a high level of quantum theory is required, which however turns out 
to be computationally too demanding for many-ion systems. Because of this, but also due 
to the weak computer power in the early days of force field development, it has become a 
common habit to parametrize the empirical force fields of ions based on single ion proper- 
ties, such as ion solvation free energies or hydration structure in small water clusters, see for 
example Ref. [l^j- Force fields based on this procedure, though, often fail to reproduce real- 
istically the ion-ion fluid structure and thermodynamics of the electrolytes at non vanishing 



concentrations, even for simple ionic solutions such as NaCl 
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14| . As has been 



recog- 
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nized in the literature, there is a strong sensitivity of solution thermodynamics 15, 
and contact ion pairing isl, to small changes in the effective pair potential between the 
interacting ions. Ion force field development is, thus, a difficult task and remains an active 
field of research 
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2l| . Recently, two studies have revisited and systematically explored 



single ion hydration free energies by scanning through a wide region of the LJ parameter 
space {cr,e} [22, l23|. This was triggered by the observation that while traditionally force 
field parameters have been adjusted in order to correctly reproduce single-ion free energies 
of solvation, one experimental observable is clearly not sufficient to unequivocally fix the two 
parameters e and a. In Hi, crystal lattice energies have been used as a second independent 



optimization target. In 



231], the single- ion solvation entropy and the effective solution ion 



size (as constructed from the ion-water radial distribution function) have been used, though 
the simultaneous optimization of two parameters, especially for the cations was problematic. 
It was seen that the three observables considered, namely the free energy of solvation, the 
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entropy of solvation, and the effective ion size, roughly matched the experimental values on 
a whole curve in the {a, e} parameter plane ,^ not allowing to single out one of the cationic 
parameter combinations as truly optimal [23]. It would be desirable to nail down the final 
cationic parameter set by benchmarking to collective, thermodynamic solution properties 
such as the electrolyte activity or osmotic pressure. Weerasinghe and Smith introduced 
and carried out this idea for NaCl solutions by reproducing Kirkwood-Buff (KB) integrals 
as determined by experiments, ensuring that a good representation of solution activity is 
obtained [2J]. The same approach has been used recently to investigate the cation-specific 



binding with protein surface charges 25|. However, the parametrization, which involves 
a free fit parameter in the mixing rule, does not conserve the free energy of solvation of 
single ions. It is therefore an interesting question, whether one can simultaneously describe 
single ion properties (such as the free energy of solvation) and ion-pairing properties by 
choosing optimized parameters for {a, e} alone or whether an additional parameter has to 
be introduced, e.g. in the form of a generalized mixing rule as was recently suggested 

An alternative method for benchmarking MD force fields has been introduced by Hess et 
al. 26| and Kalcher and Dzubiella Here, effective, MD-derived ion-ion pair potentials 

are used in a many-body corrected virial route to obtain osmotic pressures. The electrolyte 
structure at a given concentration, which forms an input to the virial equation, can be ob- 
tained directly from an MD simulation or approximately, but with much less computational 
effort, from simulations with implicit solvent l4, l26| or hypernetted-chain (HNC) integral 
equation theory [28] . It was shown that the KB derived NaCl force field 24 ] and a few alkali- 
Cl force fields proposed by Dang 29|] could reproduce experimental osmotic coeffcients in 



SPC/E water quite well, while others badly failed [26|, |27|. The reason for the failure of 
some of the force fields when considering ion-pairing properties did not become clear. 

In this work, we explore the optimization of ionic force fields based on single-ion and 
ion-pair thermodynamic properties, using the infinite-dilution solvation free energy and the 
finite-concentration electrolyte osmotic pressure as benchmarks. As for NaCl reasonable 
force fields exist, we use for Na"*" and Cl~ the parameters given by Dang 
the salts KCl, Nal, KF, CsCl, and Csl in SPC/E water and systematically vary the force fields 
of K"^, Cs^, F~, and I^. To satisfy the experimental ion solvation free energies, we confine 
our search in L J parameter space to the experimental equi-solvation free energy lines me — a 
space as calculated previously [23^. In this paper we do not vary systematically the mixing 
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30[ , and focus on 
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rule and in most simulations use the LB mixing rules. We apply the procedure proposed 
by Kalcher and Dzubiella 0] to calculate the electrolyte osmotic coefficients at a finite 
concentration of 1 M for a wide range of LJ parameters and compare to experiments. Our 
calculations are accompanied by HNC integral equation calculations that allow to efficiently 
cover and investigate a wide range of electrolyte concentrations. We systematically explore 
(a) whether ionic force field optimization is possible consistently for both single ion solvation 
energies and collective electrolyte thermodynamics without loosening the parameter space 
constraint given by the mixing rule, and (b) how the thermodynamics of electrolyte solutions 
react to a change of the LJ parameters along the equi-solvation free energy path. As a main 
result, this simultaneous optimization of the force field seems possible for the cations Cs"*" 
and , while for the anions I~ and F~ the results are less promising. This suggests 
that in the long run, and in order to consistently describe finite-concentration electrolyte 
thermodynamics, the mixing rules have to be modified and systematically optimized. By 
calculating osmotic coefficients for a solution of Csl we also investigate the transferability 
of ion parameters for Cs"*" and I~ that have been separately optimized by matching osmotic 
coefficients for CsCl and Nal. Modulo the previously mentioned restricting comments on the 
optimization of I~ parameters, ion parameters seem transferable in the sense that trends in 
the osmotic coefficients of certain ions also are found when those ions are assembled into 
different ion pairings. 



II. METHODS 



A. Simulation details 



We perform atomistic simulations using the Molecular Dynamics (MD) package GRO- 



MACS 



32| in the {N,p, T) canonical ensemble, for which the particle number A^, as well 



as the pressure p = 1 bar and temperature T = SOOA' are held constant using a Berendsen 



barostat and thermostat 



33|. The simulation box is cubic, with an edge length of L ~ 4 



nm, and periodically repeated in all three dimension and includes explicit ions and a total 
number of about 2000 SPC/E water molecules lOj. Finite size effects are not significant 
for these sizes, as shown by a previous study on similar systems and sizes 27|. The three 
dimensional particle-mesh Ewald sum is used for the electrostatics [34] with a grid spacing 
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in Fourier space of 0.12 nm in all three directions. We use an interpolation order of 4, a 
distance cutoff of 0.9 nm for the real-space interactions, and a relative strength of the elec- 
trostatic interaction at the cutoff of 10^^. Typical times for the simulations for gathering 
statistics are 150 ns for low concentrations and 40-50 ns for moderate concentrations. 

Five different salt solutions were simulated, CsCl, KCl, Nal, KF, and Csl at densities 
of 0.3 M and 1 M with 12 and 39 pairs of ions in the solution, respectively. Here, the 
concentration M denotes molarity (mol/1). For the ions, charged and non-polarizable spheres 
were used interacting with the LJ potential (Eq. [1]). For the SPC/E water the LJ parameters 
are eoo = 0.6500 (kJ/mol) and aoo = 0.3166 (nm) and are assigned only to the oxygen 
molecule of water (no parameters are related to the two hydrogen atoms). Point charges 
of qo = — 0.8476e and qh = -|-0.4238e are assigned to the oxygen and hydrogen atoms, 
respectively. For the ions we varied both parameters Eio and cxjo and present the analysis in 
Section Hill Only for sodium (Na"*") and chloride (Cl~) we used fixed parameters, those given 



by Dang 2 



30 L as they have been proven to be consistent in determining thermodynamic 
27 1 and give reasonable hydration energies 22, 23|. The LJ parameters are 



properties [25 

for Na"*", eNa,o = 0.5216 (kJ/mol) and aNa,o = 0.2876 (nm) and for CI", eci,o = 0.5216 
(kJ/mol) and aci,o = 0.3785 (nm). In this notation, the subscripts iO denote the parameters 
between ion i and the oxygen atom of the SPC/E water model. For the cross interactions 
between two ions we use the Lorentz-Berthelot mixing rules (except where noted otherwise). 

For NaCl we also tried a different force field that was optimized based on Kirkwood- 
Buff integrals 2J|. The parameters of this force field are for Na"*", eNa,o = 0.342 (kJ/mol) 
and aNa,o = 0.279 (nm) and for Cl~, eci,o = 0.547 (kJ/mol) and ac^o = 0.373 (nm). 
That force field involves a modified mixing rule for the relation between ion-ion and ion- 
water interaction, as it was not possible to fit the experimental data without breaking the 
geometric mixing rule 24 1. 



B. Effective ionic pair potentials 

We begin with a brief description of the derivation of the effective ionic (infinite dilution) 
pair potentials for the salts studied here, for more details see The pair potentials are 
derived from the radial distribution functions (rdfs) obtained within finite-concentration 
MD simulations. The rdf between a pair of atoms or ions i and j at distance r is defined 
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as gij{r; p) at a given salt concentration p. The potential of mean force (pmf) Wij{r; p) at 



35, 



36]: 



concentration p results from the rdf through a Boltzmann inversion 

f3w,j{r; p) = - ln[g,j{r; p)] = /3[<J(r; p) + w^^^ir; p)] 
where /3 = l/fc^T is the inverse thermal energy. The pmf can be decomposed into a s 



(2) 



lort- 



ranged and a long-ranged contribution, iy|[(r; p) and w'^(r;p), respectively 



37|, 



38, 



39, 



401 ] ■ The long-ranged part of the pmf is a non-specific Debye-Hiickel potential and can be 



subtracted from t(7jj(r;£) leading in this way to the short-ranged part of the pair potential 



as detailed previously 



<(r;p)=^,,.(r;p)-<^(r;p). (3) 

In the low concentration limit, the pmf between two ions reduces to their effective pair 
potential and the decomposition described in Eq. ([2]) can be written as: 

f3Vrf{r) = (3V^{r) + z,z,\b{Q)/t (4) 

where the potential is split into the short-ranged part of the pair potential, V^{r), and 
the usual Coulombic part. In this equation, the valencies for the two ion types, 

respectively, and \b{p) is the concentration dependent Bjerrum length with an infinite- 
dilution (pure water) value of Ab(0) = 0.78 nm for SPC/E water (about 10% larger than 
the real water value) 4l|. The key assumption of our derivation is now that the short- 
ranged part of the pair potential, V^J\ can be extracted from the finite concentration pmf, 
V^j'ir) ~ Wij{r]p), as calculated in (3), at not too high concentration. This is a good 
approximation, as long as the density is smaller than the density where thejiydration layers 
of ions begin to overlap, which has been found to be betwen 0.5 and 1 M [27]. On empirical 
grounds, the above procedure for the derivation of the ionic pair potentials works well for rdfs 
generated at a finite concentration of p ^ 0.3 M. Here, V^^'"(r) ~ w^J{r; p) is well fulfilled and 
accurate rdfs can be sampled with good statistics 27|. In the following, the total effective 
pair potential V^j^{r) is used as an input to the pressure calculations by the virial route and 
as an input to the HNC method. 



C. Virial route to the osmotic coefficient (l){p) 



The optimization of the ionic force fields for the salts studied in this work is done by 
comparing the derived osmotic coefficients to their experimental values. We use the virial 
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route to calculate osmotic coefficients as was done previously 27|. The osmotic pressure 
n = 2p(f){p)/i3 of the ionic solution is defined through the osmotic coefficient (f){p) at a 
concentration p. The osmotic coefficient is given through the virial equation jsol: 



9ijir;p) 



^ r ar, 

dr 



(5) 



where the indices i and j represent the two salt components and Qijir; p) needs to be evalu- 
ated at the respective concentration. 

The virial route, as implemented in this work, is not exact as it employs the infinite 
dilution pair potential V^^ . Accordingly, many-body contributions to the ion-ion interactions 
for higher densities as induced by the water are not included. It has been suggested, though, 
that many-body contributions to the pair-potential can be qualitatively included b y ta king 
into account the concentration dependence of the water dielectric constant e{p) IJ, |26 |. 



Thus, the long-range part in the pair potential V^j{r) has to be altered by using e[p) 
instead of the infinite dilution limit e(0). The following correction has been shown to lead 
to agreement of the virial route with the exact compressibility route up to a concentration 
of roughly 2 M |27|: 



[Ab(0) - \b{p)] 



(6) 



where again \b{,p) = I3e^ / ^t^SqE^p) is the Bjerrum length in the aqueous electrolyte solution 
of concentration p, and £(0) ~ 72± 1 for SPC/E water ^^-i consistent with previous studies 
4l| . The input parameter e{p) is directly calculated from the MD simulations and fitted 
through the function 



1 + Ap) 

where the values of the constant A for each salt are shown in Table [B 





CsCl 


KCl 


Nal 


KF 


Csl 


A [1/M] 


0.23 


0.24 


0.34 


0.19 


0.26 



TABLE I: Values for the A parameter in Eq. [7] as calculated from this and previous work 
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D. The hypernetted chain (HNC) approach 



A more efficient yet more approximate evaluation of the variation of the osmotic coefficient 
over a wide concentration range is possible with integral equation theory based on the HNC 
15| . The latter relates the pair correlation function gij{r) between two particles i 



closure 



and j to the pair potential in an approximative way [39|] , through: 



g,j[r) = exp 



PV:f^ + h^.ir) - c,,{r)] (8) 



where Cjj(r) is the direct correlation function, and hij{r) = gij{r) — 1 is the total correlation 
function. The HNC equation is closed by the Ornstein-Zernicke (OZ) equation of liquid 



state theory which relates hij{r) and Cij{r) 4^. For an N-component mixture with particle 
number densities p„ the OZ equation is given as hij{r) = Cij{r) + J2k=iPk I dr'^^iki^ — 
r')hjk{r'). The HNC approach uses the £(p)-dependent effective pair potential (Eq. [6]) from 
the MD simulations as input. Output is the liquid structure in form of the radial distribution 
functions gij{r). The osmotic coefficient of the salt solutions under consideration can then 

equation (Eq. [5]). This approach has been used before and details 



be calculated by the viria! 
can be found elsewhere 
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The derived osmotic coefficient-concentration curves are in 
good agreement with the MD-derived ones but start to show significant deviations above 
1.5-2M. 

We note that experiments as well as atomistic MD simulations treat the system at the 
so-called Lewis-Randall (LR) level while the implicit HNC theory uses the McMillan-Mayer 
(MM) level [sgI. At the MM level, the thermodynamic functions are calculated at constant 
chemical potential of the solvent, thus the MM and LR approaches correspond to different 
statistical ensembles. The pressure is given by the total pressure of the solution within LR 
and by the osmotic pressure of the solutes in equilibrium with the solution within MM. In 
order to compare the results between MD and the HNC approach used here, we follow the 
conversion: 

0MM = <Plr{1 + mM,)^ = ^lr"^, (9) 
where m and p are the molality and molarity of the solution, Mg the molar mass of the 



solute, po and p' the mass densities of the pure solvent and the solution, respectively [28|. 
We thus use Eq. (jO]) to convert the HNC results to the LR level. Throughout the paper the 
osmotic coefficients shown are those corresponding to the LR level, and we use the notation 
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without a subindex. 



E. Choice of parameters and optimization procedure 

As a crucial ingredient to our strategy, all Lennard- Jones parameters investigated by us 
lie on the curve that reproduces the experimental free energy of solvation for a given single 
ion, which has been calculated previously (23!]. This way, we can check whether parameter 
combinations exist that simultaneously reproduce single ion solvation as well as collective 
solution properties. All Lennard- Jones ionic parameters employed in this work are depicted 
in Fig. [TJ 

The procedure we have followed here for the optimization of ionic force fields is the fol- 
lowing: for each salt solution and each parameter set used, (a) we start with MD simulations 
of about 150 ns at a concentration of 0.3 M and (b) we simulate the same system for about 
30-50 ns at a concentration of 1 M. We determine the rdfs between the ions and extract the 
effective pair potentials from the low-concentration simulation according to the methodology 
outlined in III B[ The effective pair potential and the rdf leads to the osmotic coefficient for 
the specific solution at 1 M according to the virial route and Eq. ([5]). Note that for step 
(a) the simulation time is longer as more statistics need to be gathered for the computation 
of the pmfs at 0.3 M, compared to the rdfs at 1 M concentration. In order to obtain the 
variation of the osmotic coefficient 0(p) for a whole concentration range, we follow the HNC 
approach outlined above. 

Optimization of the LJ parameter for each ion is attempted by comparing the osmotic 
coefficient at 1 M as calculated from the MD simulations to the corresponding experimental 
values. The 0(p) curves from the HNC calculations only serve as an indication whether the 
overall behavior of is reasonably compared to the experimental curves and is not used for 
parameter optimization. 

The salts that were modeled in this work are CsCl, KCl, Nal, KF, with the goal to 
obtain optimized force fields for Cs^, K+, I~, and F^. As a check on the transferability of 
the obtained parameters, we also considered a solution of Csl with parameters optimized for 
CsCl and Nal. NaCl from Dang in SPC/E has been found to describe the osmotic properties 
and activity well when compared to experiments 



reasonable values for the solvation free energies 22 



27I . while the individual ions also yield 
23| . For this reason, we do not attempt 
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to optimize the force fields of Na"*" and Cl~ and rather consider them as a given reference. 
We begin with MD simulations of CsCl, KCl, (for which the Cl~ force field is that from 
Ref. [3o|) and Nal (for which the Na"*" force field is that from Ref. {29]). MD simulations 
are performed for those three salt solutions for all ionic parameters of Cs^, K+, and 1~ 
summarized in Fig. [HThe optimal ionic parameters for Cs^, K"*", and I~ are then estimated 
from the comparison of the resulting osmotic coefficients to the experimental values. At 
the next step, we use the optimal LJ parameters for K"*" in simulations of KF solutions and 
pursue a similar parameter space survey for F~ with the goal of finding an optimal parameter 
set for F~. The choice of the salts Nal and KF is mainly motivated by the fact that the 
standard force-fields for those salts gave very poor description of the osmotic coefficients 
in previous investigations As a consistency check, we finally take the optimized LJ 

parameters for Cs"*" and I~ and consider Csl solutions and compare the calculated osmotic 
coefficient for Csl to experiments. 



A well-known problem of ionic force fields 



431] has to be mentioned. Our MD simulations 



show that very low Sio parameters for the cations (Cs+ and K"*") lead to unphysical clustering 
of the ions for CsCl and KCl even at low concentrations of 0.3 M. Accordingly, no pmfs 
could be extracted and the corresponding parameters (Csl, Cs2, Kl, K2, K7, and K9) are 
neglected for the calculation of the osmotic coefficient and for the optimization procedure 
of the Cs"*" and K"*" force fields. On the other hand, very low Eio for the anions (1~ and F~) 
do not lead to similar clustering for Nal and KF, respectively, and can still be considered 
as candidates for an optimized force field. Interestingly, ion clustering was observed for all 
salts at a concentration IM in some regions of the LJ parameter space. Specifically, the 
ions in KF clustered for the parameter sets F12, F13, for which the Kll parameter was 
used. For Fl, F2, we could not obtain reasonable results, as the ion- water system could 
not be energetically optimized within MD, thus those parameters had to be eliminated as 
well. Apart from these restrictions, we have used all other parameters in Fig. [1] for the MD 
simulations. Parameter sets for which results are not shown for the osmotic coefficients in 
Fig. [5] are the ones that led to clustering of the ions in the aqueous environment. 
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III. RESULTS AND DISCUSSION 



A. Structural properties 

1. Radial distribution functions 

A robust measure of electrolyte structural properties is the ion-ion radial distribution 
function (rdf), which shows distinct structural signatures that differ among different salt 
solutions. Here, we have calculated the rdfs for all systems at two concentrations, 0.3 M 
and 1 M, respectively. The rdfs for the two concentrations and the same LJ parameters do 
not differ significantly, apart from the height of the rdf peaks. The heights of the contact 
and the solvent-separated peak indicate different hydration properties of the ions. All curves 
exhibit strong electrostatic screening and reach the asymptotic value of 1 below a distance of 
2 nm for the 1 M concentration. This is consistent with the small screening length of about 
0.25 nm at 1 M. As shown in Fig. O for CsCl and KCl the first peak in the cation-anion 
rdfs is much higher than the second one, indicating close contact of the anion and cation 
and predominant direct ion pairing. For KF and Nal, the first peak in the cation-anion 
rdfs is of similar height or lower than the second one, indicating that water enters between 
the anion and cation indicative of indirect ion pairing. Due to the variety of LJ parameters 



used in this study, the re^ 
salts shows rich behavior 



ative strength of direct and indirect ion pairing for the different 
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19| . The height of the first rdf peak at 0.3 M for CsCl ranges 
from about 6 to 30, for KCl from 5 to 25, brought about by variations of the cationic force 
field parameters, as illustrated in Fig. [21 For the ion pairs that show pronounced solvent- 
separated pairing, Nal and KF, the influence is much less, here the height of the first rdf 
peak varies for Nal from 2.5 to 3.5 and for KF from 9 to 12. Note that these variations are 
induced by changing the anionic force fields. This demonstrates that by changing force field 
parameters that keep the single-ion properties invariant (as judged by the ion-water rdf or 
the solvation free energy), the ion-pairing properties can be significantly affected [l3]- The 
position of the contact peak also shows interesting behavior. As an example, we present the 
trends for some representative LJ parameters for the KCl cation-anion rdfs: for the order of 
the bare LJ radius axoi (^ko{.K11) < aKo{K8) < aKo{K6) < aKo{K5) < aKo{K3), the 
ordering of the contact peak position between K"*" and Cl~, t'^^qi is r]^^i{K3) < r'^(ji{K5) < 
''^KCii^^) ^ ''^Kcii^^) ^ ''^Kcii^^^)^ visible from Fig. O We find similar trends for 
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CsCl, Nal, and KF, as well. Ions with smaller bare Lennard Jones radius thus show larger 
ion-ion separation, which clearly has to do with the fact that the interaction strengths {sio) 
of the single ions are varied along with cTjo- In addition to the peaks of the rdfs we also study 
the minima in the rdfs. The positions of the first and second minima in the rdfs, ri and r2, 
and the distance at which the rdfs vanish, rg, are given for representative salt parameters 
at a concentration of 0.3 M in Table HTl 





Cs6Cl 


Cs9Cl 


K5C1 


KUCl 


Nail 


NaI4 


K11F5 


ro (nm) 


0.29 


0.30 


0.27 


0.28 


0.28 


0.28 


0.29 


ri (nm) 


0.43 


0.43 


0.40 


0.40 


0.41 


0.40 


0.43 


12 (nm) 


0.67 


0.68 


0.63 


0.63 


0.64 


0.65 


0.66 



TABLE II: The rdf characteristics at 0.3 M for the optimized ion parameters; see Tab. III. The 
distance (ro) at which the rdf vanishes to zero and that of the first (ri) and second (r2) minimum 
in the cation-anion rdfs for representative parameters of Cs"*", K^, I~, and in CsCl, KCl, Nal, 
and KF aqueous environments, respectively. 



We have also compared the g{r)s from our MD simulations with results from HNC for a 
few different salt parameters and observed only small differences; see Fig. |3l The reason for 
the deviations is the approximate treatment of statistical mechanics in HNC. The differences 
are mostly visible in the height of the first and second peaks in the g{r). 



2. Short-ranged potentials of mean force 

We derive short-ranged pair potentials for all salt parameters that do not lead to ion 
clustering from Eq. [3l Examples are shown in Fig. |H In accordance with the rdfs, the 
pair potentials for Nal and KF reveal deeper second minima, while the first two minima for 
CsCl and KCl show roughly similar depth. The second minimum of KF is broader than for 
the other salts, indicating an unusual water configuration in the solvent-separated ion pair. 
The Cs^-Cs^ and K"''-K"'" pair potentials show smaller oscillations compared to the Na"*"- 
Na"*" potential. This has been observed before, and was rationalized by the fact that Na"*" 
has tightl y b ound hydration shells giving rise to energy barriers when two sodium cations 
approach 27| . In the case of the anion-anion pmfs, the F~-F~ pmf shows deeper minima and 
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stronger oscillations than the C1~-C1~ or the potential. Going from the small fluoride 
to the large iodide, one observes a trend towards a soft-sphere like potential. Similar to 
the rdfs, we do not observe a simple dependence of the position and depth of the minima 
in the pmfs with the variation of the LJ parameters, again due to the simultaneous change 
of both the LJ radius and interaction strength along the lines of constant solvation free 
energies. Note again that there is substantial variation in the cation-anion potentials for the 
different force-field parameters, which gives hope to be able to fit the osmotic coefficients. 
The scattering in the Cl^— Cl^ potentials (note that the Cl~ force field is not changed) is 
due to bad sampling statistics as the ions typically do not get very close. 



B. Osmotic coefficients 



Having determined the rdfs and short-ranged pair potentials for different LJ parameters, 
we next calculate the osmotic coefficient using the virial route based on the rdfs and 
pair potentials from MD simulations (not HNC) as explained in Section [Tll The results 
for CsCl, KCl, Nal, and KF at a concentration of 1 M are summarized in Fig. [5] where 
we have added lines as guides to the eye in order to bring out the main features of the 
results more clearly. The error bars for the calculated are in the range ±0.01-0.05. The 
experimental values of the osmotic coefficients 4^, |45|] for each salt are also shown, on which 
we base our optimization of the ionic force fields. Inspection of this figure reveals that for 
the salts CsCl, KCl, and KF shows a maximum for intermediate values of Eio- For Nal, 
no significant variation of with Sio is observed. We first focus on the cations Cs"^ and K"*" 
in CsCl and KCl, respectively, and panel (a) in Fig. [51 As more clearly shown by the lines 
added as guides to the eye for the CsCl and KCl data, there are two crossings between the 
simulated and the experimental values for 0, so in principle there are two optimal force field 
sets for these cations. Note that we also included for the Dang parameters for KCl and 



CsCl 



30j as open symbols. Interestingly, though they are somewhat off from the curve 



corresponding to the experimental solvation free energy (see Fig. [T]), they show quite good 
agreement with the experimental data for within the error. 

Inspection of panel (b) reveals that the situation for the anions I^ and F~ is distinctively 
different. For all iodide parameters, the corresponding values are close to the experimental 
value (within the error), but show very little variation. This suggests that varying the force 
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field for I lias no considerable effect on the osmotic coefficients. A very similar result was 

n 

obtained for chloride in NaCl in previous simulations [24]. However, this does not seem to 
be generally true for anions, as KF in the same panel shows a much larger variation in 0, 
ranging from about 0.4 up to 0.7. For most of the results for KF shown in Fig. MJo), the 
optimized Kll parameter set from Table UTTl was used. However, as seen in the figure, the 
use of the equally optimal force field K5 does not lead to qualitative changes. For all F~ 
parameters, the osmotic coefficients of KF are too low compared to the experimental osmotic 
coefficient and even for the best F~ parameter combination are too low by about 0.25, which 
is larger than the error bars. We conclude that while force fields for Cs"*" and K"^ can be 
derived with relative ease, the situation is more complicated for the anions considered by 
us: while works quite well (which seems to be coincidental since the variation of (p with 
Eio is very small), the optimization for F~ fails though here the variation of (p with ejo is 
pronounced. 

Using the HNC approach, we have also calculated the variation of the osmotic coefficient 
with the concentration for all salts and LJ parameters investigated here. Curves for some of 
the parameter sets used are shown in Fig. [6]together with experimental data. For some of the 
LJ parameters for fluoride, (p for KF diverges above a certain concentration, for those cases 
the corresponding curves are cut above the concentration of 2M. Inspection of the overall 
shape of the curves reveals that there is good agreement with the experimental curves for 
some LJ parameters for CsCl, KCl and Nal, in contrast to KF. This is not unexpected, 
since for KF simulation results for did not match the experimental value at 1 M for any 
parameter combination. Note that there are small differences between the HNC results 
(lines) and MD results at 0.3M and Im (symbols) for 0, which are of the order of the error 
inflicted by the virial approximation and the simulation numerical scatter of about ± 0.05, 
see our previous discussion |28| . 



C. Optimum ionic force field parameters 



Based on the MD osmotic coefficient results, we suggest optimal LJ parameters for Cs"^, 
K^ for use with Cl^ parameters from Dang 30|] , and optimal parameters for I~ for use with 
Na"*" from Dang 29|; we also suggest a force field for F~ derived from simulations of KF 
using our optimized force field for K"*". For most ions, two parameter sets are suggested and 
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summarized in Table IIIII For the cations, Cs"*" and K"*", we take two values closest to the 
crossing points of the fitting curves for (j) with the experimental data, see Fig. [5], namely Cs6, 
Cs9, K5, and Kll. For K"*" a slight ambiguity exists, as for example also the K8 force field 
matches the experimental data. We rather chose Kll because that parameter is further 
away from the parameter sets K7 and K9 for which the ions were found to cluster. For 
iodide, our choice of parameters is less obvious; basically all parameters within the error 
bars are equally acceptable. 



ion /label 


Gio (nm) 


ejo(kJ/mol) 


Cs6 


0.333 


0.5 


Cs9 


0.325 


1.0 


K5 


0.31 


0.41 


Kll 


0.293 


1.26 


11 


0.45 


0.1 


14 


0.425 


0.32 


F5 


0.3665 


0.1 



TABLE III: The optimized Lennard- Jones parameters for Cs"*", K^, I , and F as extracted by 
the MD simulations and the comparison of the osmotic coefficients with experimental data. 

For KF, we performed two distinct sets of MD simulations, the first with the Kll force 
field and the second with the K5 parameters. Interestingly, both K+ give comparable results 
for 0, see Fig. [5)3, meaning that the redundancy found with KCl seems to be also present 
for KF. However, none of the KF parameters reproduces experimental values for 0, so we 
chose a single force field for F" that shows the least deviation. 

In Fig. m pmfs for the optimal parameters are compared with a pmf of a non-optimal LJ 
parameter set. As expected, for Cs^ and K+ the cation- anion pmf of the non-optimal force 
field shows larger deviations from the optimal force field results; for I~ all anion-cation pmfs 
are quite similar. 

We briefly return to the discussion on peak heights and ion pairing [isl, [3]. For the 
optimal force fields in Table IIIH we find for the height of the first contact peak in the rdfs 
at 0.3 M concentration, values of about 9.5 and 5.5 for K5C1, and KllCl, respectively; 14.8 
and 8.6 for Cs6Cl and Cs9Cl, respectively; about 2.3 for both Nail and NaI4, respectively, 

16 



and 3.3 for K11F5. The order of the peak height is KCl ~ CsCl > KF > Nal, consistent 
with previous theories on ion pairing 13], according to which the tendency to form direct 
ion pairs goes down as the ion sizes become more dissimilar. We note that this ordering of 
contact pair formation probabihty is only realized for the optimized force fields, non optimal 
force field combinations can easily lead to partial or complete reversal of this ordering. 



D. Transferability of the optimized ionic force fields 

We so far were occupied with finding force field parameters that match experimental os- 
motic coefficients best. We now turn to a separate question and check whether the optimized 
force fields presented in the previous section are transferable. To that end, we perform a set 
of MD simulations for Csl in water, for which the parameters for Cs^ and are taken to 
be the optimized force fields from Table UTTl The hope can be not be to match experimental 
osmotic coefficients for Csl perfectly, as the Iodide parameters are not perfect by themselves 
(when compared with experimental osmotic coefficient data for NaT). Rather, we intend to 
check whether the trends of the simulated osmotic coefficients of Csl reflect the properties 
of the Cs"^ and I~ force fields. The specific salts modeled are Cs6I4 and Cs9I4 (we did not 
check the II parameter set, as variations of iodide parameters do not seem to considerably 
affect the osmotic coefficients, as was seen from Fig. Ofor Nal). The MD and HNC results 
for the osmotic coefficient are shown in Fig. W^a). Note that for Cs6I4 we only show the MD 
data for 0.3M, since the simulation data at IM show bad convergence behavior, probably 
due to the vicinity to the crystalization transition for this particular force field (note that ex- 
perimentally Csl has the smallest maximal solubility of all considered salts in this paper, of 
about 3M, so such problems might be anticipated). We find overall good agreement between 
the MD and HNC results with deviations in smaller than 0.05. There are quite sizeable 
deviations between the theoretical predictions and the experimental data. However, note 
that we have used the optimized parameters for both Cs"*" and I~, based on osmotic coeffi- 
cients for CsCl and Nal, and that the deviations from experimental data for Csl in Fig. [Tl^a) 
are of the same order as the deviations observed for CsCl and Nal in Fig. [51 Furthermore, 
the deviations from experimental data go in the expected direction, namely, the simulation 
prediction for Cs6I4 lies below the experimental value, while Cs9I4 lies above (similar to 
the data for Cs6Cl and Cs9Cl in Fig. [5ti). We tentatively conclude from these data that 
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the force fields obtained by our optimization strategy are transferable, meaning that if one 
fixes the force fields of ions and by optimizing the osmotic coefficients of the salts 
AB and CD, the osmotic coefficient of the salt AD should come out approximately right 
without further adjustment. Such reasoning of course assumes that one can optimize the 
salts AB and CD and therefore does not apply to F", where the optimization of KF fails in 
the first place. So one sees that while optimization and transferability are logically distinct 
operations, transferability presumes successful optimization of force field parameters. We 
will come back to this point in the conclusions. 

As an additional check of our methodology, we also used a different force field for NaCl 
that was designed to reproduce experimentally determined Kirkwood-Buff integrals (and 
thus experimental activity coefficients) 2J]. We performed simulations with those NaCl 
parameters using the mixing rules proposed, namely geometric mixing including a freely 
adjusted mixing coefficient (and not the Lorentz-Berthelot rule which we were using up to 
this point). The MD and HNC results for the osmotic coefficient are shown in Fig. [7](b), 
together with the experimental data and those derived previously from the NaCl Dang pa- 



rameters 
fields 



27, 



29 



301]. It is evident from this figure that the Kirkwood-Buff NaCl force 



2J] give good agreement with experimental values for the osmotic coefficient, indi- 



cating that activity coefficient and osmotic coefficient data probe the same ionic features, 
namely the pairing characteristics in aqueous solution. 



IV. SUMMARY AND CONCLUSIONS 



We have used a combination of MD simulations and statistical mechanics analysis to 
systematically explore and optimize force field parameters for ions in aqueous solution, with 
particular emphasis on the interplay of single-ion and ion-pair thermodynamics. 

The ion parameters considered by us, namely the Lennard- Jones radius and strength, 
were confined to a curve on which the experimental single ion solvation free energy is re- 
produced {23! • For a whole number of specific force fields on those curves, we constructed 
effective ionic pair potentials which led, through the virial route, to electroljd;e osmotic co- 
efficients. These were then compared to experimental values at 1 M concentration. We 
have used the MD-virial route derived osmotic coefficients to optimize the LJ parameters 
for the single ions Cs"*", K"*", I~, and F~ based on data for CsCl, KCl, Nal, and KF, respec- 
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tively, treating Na"*" and Cl~ as reference ions. This optimization strategy works well for 
the cations Cs"*" and K"*", for which we obtain, due to the peculiar shape of the MD-based 
osmotic coefficient curves, two optimized force fields for each ion. This proves that at least 
for the cations, the simultaneous description of single-ion and ion-pairing thermodynamics 
is possible if one systematically explores the full LJ parameter space without the need to 
modify the combination rules. For anions, on the other hand, the optimization is more 
problematic. We could not get reasonable match of experimental osmotic coefficient data 
varying the LJ parameters of F~ in KF, simply because the maximum in the simulated 
osmotic coefficient is significantly below the experimental value. Iodide can be more or less 
satisfactorly optimized based on Nal data, but this seems to be mere coincidence, since 
NaT osmotic coefficients shows very little dependence on the LJ parameter, which indicates 
a basic limitation of our approach. With these restrictions in mind, we checked for the 
transferability of our optimized force fields by considering the osmotic coefficient of a Csl 
solution, i.e. an ion combination that was not targeted in the optimization process. Based 
on the far-from-perfect performance of the I~ force field in NaT, we judge the transferability 
properties of the force fields as satisfactory since in particular the trends of the CsCl osmotic 
coefficients upon variations of the Cs^ force field are fully recovered when looking at the 
trends of the Csl osmotic coefficients. 

These results suggest that for truly accurate non-polarizable ion parameters, one appar- 
ently needs to lift the constraint of the simple mixing rules and introduce an additional 
scaling parameter in the mixing rule, in agreement with previous approaches 2J, |25|]. How- 
ever, it seems desirable and possible to introduce such a mixing parameter only in the 
cation-anion interaction, so that the anion-water and the cation-water interactions stay the 
same as used in the single-ion solvation free energy optimization (the cation-cation and 
anion-anion interactions are much less important). That way, one would correctly describe 
the single-ion properties (as embodied in the infinite-dilution solvation free energy) as well 
as ion-pairing properties (as important for osmotic and activity coefficients). As an un- 
fortunate byproduct, one would have to optimize this mixing parameter for each ion pair, 
requiring substantial simulation efforts. 
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FIG. 1: All Lennard-Jones Eio and ctjo parameters for the ions (a) Cs"*", (b) K"*", (c) I , and (d) 
used in the optimization procedure in this work. All parameters lie on the curves on which 
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FIG. 2: Cation-anion radial distribution functions (rdfs) at a concentration of 0.3 M for various 
parameter sets for (a) CsCl, (b) KCl, (c) Nal, and (d) KF. The Cl~ and Na^ force fields are taken 



from Dang 



301], the K+ force field in (d) is Kll. 



24 




I I I I I I I I I I I I I I I I I I I 

0.2 0.4 0.6 0.8 1 1.2 

r (nm) 

FIG. 3: Radial distribution functions (rdfs) at a concentration of 0.3 M for the Cs6 parameter set. 
Red (solid) lines correspond to the MD results and blue (dashed) lines to the rdfs as generated 
through the HNC approach (based on the pmfs from MD). From top to the bottom, the rdfs for 
the cation-anion, cation-cation, and anion-anion are shown. 
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FIG. 4: Short-ranged potentials of mean force (pmfs) at a concentration of 0.3 M for the optimal 
LJ parameters (see Section llll Cp for (a) Cs (sets Cs6, Cs9), (b) K (sets K5, Kll), (c) I (sets II, 
14) and (d) KF (set K11F5). In each panel, from top to the bottom, the pmfs for the cation-anion, 
cation-cation, and anion-anion are shown. For comparison, in all panels the pmfs for a non-optimal 
parameter set is included: (a) set Cs3, (b) K3, (c) 17, and (d) K11F3. 
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FIG. 5: The osmotic coefficient (p as calculated from the MD simulations by means of the virial route 
for CsCl, KCl, Nal and KF in a 1 M solution. In (a) the magenta(dashed) and turquoise (dotted) 
lines represent the experimental results for CsCl and KCl, respectively. In (b) the green (dotted) 
and orange(dashed) lines are the experimental values for Nal and KF, respectively. Representative 
results for KF for the K5 parameter set are also shown. The open symbols are the osmotic 
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FIG. 6: The osmotic coefficient (\) as calculated through the HNC approach (lines, see text) of 
representative LJ parameter sets for (a) CsCL (b) KCl, (c) Nal, and (d) KF as a function of 
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45l | are also shown. The symbols are the values 



directly from MD for the parameters sets proposed in Table IIIII In panel (d) all KF results are 
obtained by using the Kll parameter set for the potassium ion. 
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FIG. 7: MD results (symbols) for the osmotic coefficient, plotted together with the HNC results 
(broken lines) as a function of the salt concentration (in M). Panel (a) shows results for Csl 
and serves as a check on the transferability of the proposed force fields for Cs"*" and I~. Panel 
(b) depicts the results for NaCl from the Kirkwood-Buff derived force field (kbff) [24] and the 
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